ANALYSIS AND COMPUTATION OF A DISCRETE 
KDV-BURGERS TYPE EQUATION WITH FAST DISPERSION 
AND SLOW DIFFUSION 



ZVI ARTSTEIN, C. WILLIAM GEAR, lOANNIS G. KEVREKIDIS, MARSHALL SLEMROD, 

AND EDRISS S. TITI 



Abstract. The long time behavior of the dynamics of a fast-slow system 
of ordinary differential equations is examined. The system is derived from 
a spatial discretization of a Korteweg-de Vries-Burgers type equation, with 
fast dispersion and slow diffusion. The discretization is based on a model 
developed by Goodman and Lax, that is composed of a fast system drifted by 
a slow forcing term. A natural split to fast and slow state variables is, however, 
not available. Our approach views the limit behavior as an invariant measure 
of the fast motion drifted by the slow component, where the known constants 
of motion of the fast system are employed as slowly evolving observables; 
averaging equations for the latter lead to computation of characteristic features 
of the motion. Such computations are presented in the paper. 
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1. Introduction 

We consider the slow-fast system of ordinary differential equations 

, . dUk Uk{Uk+i - Uk-i) _ Uk+i - 2Uk + Uk-i 

^ dt ^ 2h ^ 

k — 1, 2, • • • , iV, with periodic boundary conditions Uk = Uk+N, and where h > 
and e > are small parameters, such that e <C /i in a sense that will be specified 
later. It is an elementary observation that p.ip is a spatial discrete analogue of the 
classical Burgers' equation 

(1.2) ut + uux = euxx, 

subject to periodic boundary conditions, with a basic domain [0, 27r], and h — ^ 
is the size of the uniform mesh with N points. But a more careful analysis (see, 
for example, Goodman and Lax [11]) shows that a higher order approximation of 
(|l.ip is provided by the Korteweg-de Vries-Burgers' type equation 

(1.3) Ut + u{Ux + —Uxxx) = £Uxx: < X < 27r, 
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with periodic boundary conditions 

(1.4) u(0,t) = u(27r,t). 

In section 2 we establish a formal derivation of p.ip from (|1.3p , following Goodman 
and Lax [11]. This relation of (|l.ip to (|1.3p . under the assumption that e h, 
reflects a competition between fast dispersion and slow diffusion mechanisms. Our 
goal in this paper is to find an effective equation governing the long time behavior 
of the "slow" diffusion, obtained via an "averaging" of the "slow" right hand side 
of (II. 1|) against the fast oscillations originating from the left hand side of (jl.ip . 
Intuition suggests that averaging will be relevant; however, implementing classical 
averaging results is not possible; we elaborate on this in the body of the paper, 
but note here that the source of the difficulty is to establish the existence and the 
identity of the "angle variables" in, say, a Hamiltonian representation, with respect 
to which averaging should be carried out. Solutions of (jl.ip can be presented as 
images of solutions of a Hamiltonian system, but the latter are not bounded, so do 
not induce relevant information (this issue is discussed in section 3). Our approach 
relies on an alternative new theory of computing slow observables for singularly 
perturbed differential equations [3]. This theory has the distinct advantage of not 
requiring an explicit knowledge of the "angle variables" and for this problem will 
produce a system of equations for the slow observables of our approximate system 
of ordinary differential equations. The Lax invariants, which are traces of powers 
of a matrix component in a Lax pair, as developed systematically in Goodman and 
Lax [11] and the related work of Kac and van Moerbeke [13], are slow observables. 
We employ these and other invariants of the fast flow in the computations of the 
solutions of the discrete approximation. 

The implementation of averaging techniques for systems of equations which are 
separated into fast and slow components is well known. These techniques in addi- 
tion to component splitting often require explicit information on the fast dynamics, 
e.g. time periodicity, existence of limit cycles, or known stochastic behavior (see for 
example the survey of Givon, Kupferman, and Stuart [10], see also E et al. [8] and 
Kevrekidis et al. [14]). The issue of averaging in the absence of explicit knowledge 
of the fast dynamics was first considered by Artstein and Vigodner [7] in 1996 and 
continued in Artstein [1] , [2] , Artstein and Slemrod [5] , [6] . A rigorous theory of av- 
eraging in the absence of both a component split and explicit limit characteristics of 
the fast dynamics was first given in [3] . The possible application of theory of [3] to 
a system with fast dispersion and slow diffusion was alluded to in [3] and this paper 
completes part of that program. In fact when the theory is coupled with techniques 
such as projective integration it provides, as we demonstrate in section 6, a fast 
and accurate method for computation of fast dispersion systems undergoing a slow 
perturbation. Of course our method does require some special features of the fast 
system, i.e. a rich family of first integrals of motion (invariants) of the type usually 
associated with, but not limited to, for instance, completely integrable Hamilton- 
ian systems. System (|l.ip is particularly interesting in that as far as available in 
the relevant literature, e.g. Moser and Zehnder [17], there is no known one-to-one 
correspondence between its fast part (when e = 0) and a Hamiltonian system, nor 
do we know whether solutions of this fast system can be presented as the image of 
solutions of a Hamiltonian system with bounded dynamics. This fact renders the 
splitting into fast and slow components difficult at best and perhaps impossible. 
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The paper has five sections after this Introduction. Section 2 derives the semi- 
discrete approximation to the system in a manner suggested in the paper of Good- 
man and Lax [11]. It aUows us to represent, in Section 3, the fast part of the system 
in terms of a Lax pair of matrices, thus producing a set of invariants of the fast dy- 
namics. We also comment there on the Hamiltonian-type structure of the system. 
In section 4 we show that the positive orthant is invariant under the dynamics of 
the full system (jl.H) . Moreover, we use this fact to show the global existence of 
solutions to (jl.ip with initial data in the positive orthant. In Section 5 we apply 
the theory of [3], displaying Young measure solutions to the limit problem, and the 
resulting evolution equations for the slow observables of our equation. Finally, in 
Section 6, we present numerical results for the relevant "effective slow" dynamics 
obtained in the limit as the parameter e —^ 0. 

2. Derivation of the semi-discrete system 

In this section we establish a formal derivation of (II. ip from ()1.3p using the 
dispersive approximation scheme of Goodman and Lax [11] suitably modified for 
our viscous problem; for convenience of comparison we use the notation in [11]. 

For integers k set Uk{t) = u{kh,t). Hence periodicity of m in x implies Uk is 
periodic in k with period N where N = Assume h is chosen so that N is an 
even integer. Taylor's theorem yields 

(2.1) Uk+i - C/fc = hua,{kh, t) + ^h^u^^ikh, t) + ^h^u^xx{kh, t) + 0{h'^) 
and 

(2.2) Uk - Uk-i = hu^{kh,t) - lh^u^.^{kh,t) + \h^Urx^{kh,t)+0{h'^). 
Addition of ((2T|) and ((2?2|) yields 

(2.3) Uk+i - Uk-i - 2hu^{kh, t) + Ih'^u.^.ikh, t) + 0(/i4), 



while subtraction of (|2.2p from (|2.3p gives 

(2.4) C/fc+i - 2Uk + Uk-i = h\,,ikh, t) + 0{h^). 

Notice that the right hand sides of the latter two equations correspond to the 
expressions which constitute the continuous equation (|2.ip . Thus the continuous 
system (|1.3p is formally Oih?) equivalent to the semi-discrete, i.e., continuous in 
time only, discrete and iV-periodic in fc, system 

dUk , Uk{Uk+i-Uk-i) Uk+i-2Uk + Uk-i 

^ + 2h =' ' 

with initial condition 

(2.6) UkiO) =uo{kh), k^l,...,N 

(which is our original system p.ip ). The A^-periodicity means that in (2.5), and 
in equations later on, we use Ug = Un and Un+i = Ui. It is important to notice, 
though, that here the step size of the approximation and the dispersion coefficient 
are changed in a correlated way through the equality hN = 2tt; thus, a finer grid, 
namely, a larger but fixed A^, corresponds to a smaller dispersive coefficient; yet we 
always take e h. 

Goodman and Lax [11] examined the limit as /i — > of the non- viscous version 
of (jl.Sp via the analogous semi-discrete approximation. They proved, in particular. 
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that during the time period prior to the formation of a shock in the KdV type equa- 
tion, the solutions of the discrete approximation converge uniformly to the solution 
of the continuous one; but after the time instant when a shock forms in the Burgers 
equation, the discrete approximation generates oscillations. We consider equations 
(|2.5p - (|2.6p as the semi-discrete approximation of the KdV-Burgers equation (jl.Sp : 
we wish to study the latter for h fixed and as e — > 0, namely, on long time intervals. 
To this end we rewrite (|2.6p in the time scale r = et, yielding 



(2.7) 



dUk UkiUk+i - Uk-i) _ Uk+i - 2Uk + C/fe-i 
dr 



2he /i2 
Denote by U the vector {Ui, - ■ ■ , Un)', system (I2.7P can be rewritten as: 

dU 1 



(2.8) 
where 

F(f/) = 



dr e 



FiU) + GiU), 



( Ui{U2 - Un) \ 

U2{U3-Ui) 



2h 



( U2~ 2Ui + Un \ 
Us - 2U2 + Ui 



and G{U) = 



1 

7? 



\U,- 2U^ 



N 



u. 



N-1 



\ Un{U, - Un-i) j 

Now, in (|2.8p (equivalently in (|2.7[) ) we identify two contributions: The "fast part", 
namely the vector field e^^F{U) (for h fixed), and the slow part, namely, G{U), 
which corresponds to the diffusion. 



3. Analysis of the fast equation 

In this section we examine the "fast part" of the vector field determined by (|2.7p . 
i.e. we consider (I2.7P after replacing its right hand side by zero. We seek to find 
invariants of the motion. These are not affected when multiplying the equation by 
a constant. Hence we write now the fast part of equation (|2.7[) suppressing the 
coefficient (2ft,e)~^; this amounts to a change of time scale, a — {2he)~^T, without 
affecting other features of the dynamics; namely, we examine the equation 

(3.1) ^ + Uk{Uk+i-Uk-i)^0, 

where Uk{0) > and k — 1,2, . . . , N with N even and with periodic boundary 
conditions (in particular, in (|3.ip we interpret Un+i — Ui and Uq = Un)- 

Observation 3.1.T/ie product U1U2 ■ ■ ■ Un is an invariant of equation \8.1]) : 
also the product C/2C/4 • • ■ Un of the even coordinates and the product U1U3 ■ ■ ■ Un-i 
of the odd coordinates of the state, are invariants of equation \3.1\) . (We, actually, 
arrived at the above quantities guided by numerical simulations). 

Proof. Straightforward, by showing that the time derivative of each of the terms 
is equal to zero (in fact, establishing invariance for two of the quantities implies 
that the third is also invariant). 

Corollary 3.2. The positive orthant, determined by the condition Uk > for 
all k, is invariant under the dynamics of (3.1) (later in the section we establish a 
stronger property). 
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In the closing paragraph of this section wc identify Hamihonian-type properties 
of the fast system p.ip . Now we identify ^ + 1 first integrals of the fast dynamics, 
namely, y + 1 constants of motion. To this end we review the results of Goodman 
and Lax [11] for the fast system, equation (|3.ip . and also borrow extensively from 
the presentation of Kac and van Moerbeke [13] and Moser [15], [16]; as far as the 
authors can tell the first representation of (13. 1|) in terms of Lax pairs was given 
by Moser [15], but, as was noted by Moser himself, it was analogous to that of 
Flaschka [9]. Moser, however, imposes Dirichlet boundary conditions while here 
we are interested in periodic boundary conditions; this case was also discussed in 
Moser and Zehnder [17]. 

We make a change of variables 

(3.2) Al = Uk. 

With the additional change of time scale s — the system (|3.ip is transformed 
into 

(3.3) ^^A,iAl^,~Al^,), k^l,...,N, 

with periodic boundary conditions, hence in (I3.3|l Ai — Ajv+i and Aq ~ A]\j. We 
note that any positive solution of the system (|3.ip can be rewritten as a solution of 
p.3p : by Corollary 3.2 the product A1A2 . . . An is also an invariant of the equation. 
Hence the positive orthant in is invariant with respect to (j3.3p . and equations 
p.ip and (|3.3p are equivalent for solutions with positive initial data. 

Now set A — {Al, . . . , Aj^), an iV-dimensional vector of real numbers, and extend 
Ak periodically by Ak+N = ^fc- 

Define the operation T on elements A as the translation {TA)k = Ak+i and 
use the abbreviation TA = A+; likewise, is the translation in the opposite 
direction and T~^A = 

Following Goodman and Lax [11], Kac and van Moerbeke [13] and Moser [15], 
[16], we associate with each vector A the operators 

(3.4) L = AT + A_T-^, 

where A acts as a multiplication operator, i.e., when T is written in its matrix repre- 
sentation then AT means that the i-th row of T is multiplied by the corresponding 
element, namely Ai, of A. The matrix L is symmetric with diagonal elements equal 
to 0. For the convenience of the reader and further reference we give below the 
explicit terms when = 6; it should make the general case transparent. 
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With the aid of L we now find conserved quantities for the equation (|3.ip . To this 
end define 



(3.6) 



B = AA+T^ - A-A.-T^^ , 
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which is anti-symmetric. Here A = T ^A- = T ^A. Furthermore, the commu- 
tator [B,L\ = BL — LB is given by 

(3.7) [B, L] = A{Al - Al)T -A-{A- A—)T'^. 

Note that (|3.3p is equivalent to 



(3.8) ^ = [B,L] 



Then, of course. 



ds 

dL^ dL dL 



and more generally, 

(3.10) ^ = IB,LP] 

ds 

for any natural number p. Since the right hand side of p.lOp is a commutator its 
trace is equal to zero. Therefore, for each p we have 

dLP 

(3.11) tr——=0 

ds 

where tr signifies the trace. For p odd the trace of LP is zero, so it does not reveal 
any information concerning equation ()3.3p . But for p even, p = 2,4, . . . ,iV, the 
term trL^ is a nontrivial polynomial of degree p. Thus, we have arrived at non- 
trivial conserved quantities, namely, first integrals, of equation p.3p . For example, 
the first two non-trivial traces for N > 6 are 



N 
fe=l 



N 



(3.12) trL^^2j2Al 
and 

(3.13) trL^ = E(24^'+i + + ^ti)')- 

fc=i 

The last trace for the case = 6 is 

6 

(3.14) trL^ = ^(A2(Ati+4 + A2+i)2 + 4_^A2A2+i) + 12AiA2A3A4A5A6. 

fc=i 

Remark 3.3. As stated in Observation 3.1 the product of the Uk is conserved 
by the equation iS. 1]) : by i3.13\) each of the stays uniformly bounded. Hence, for 
an initial condition with positive coordinates, the quantities Uk along a trajectory 
are bounded and bounded away from zero; in particular, the solution exists for all 
time. 



Remark 3.4. Invoking the relation \3.2{l . the ^ invariants for i3. 3\) give ^ 
invariants for i3.1\) . It is not easy to verify analytically that the ^ invariants are 
functionally independent, namely, that they identify distinct ^ dynamics; numerical 
simulations show, indeed, that their gradients are linearly independent, indicating 
their functional independence. In addition, and as noted in Observation 3.1, the 
following quantity is also an invariant of i3.1\) . 

(3.15) U1U2 ---Un. 
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It has been shown by simple algebra for N = 6 that the other invariants mentioned 
in Observation 3.1, namely 

(3.16) U2U4---Un and C/1C/3 • • • C/w-i 

are not independent of the ^ + 1 noted above and that this is hypothesized to be 
true for all even N . Altogether we get ^ + i functionally independent invariants. 
We do not know whether there are any independent global (i.e. not restricted to 
a specific lower dimensional manifold) first integrals in addition to the + 1 ^^^^ 
we have identified; as we explain later on, the numerics indicates that the list is 
exhausted. 

Considerations of the number of invariants may be significant as, under an appro- 
priate transformation, solutions of (|3.ip are differences of solutions of a completely 
integrable Hamiltonian system. This was noted in Goodman and Lax [11], refer- 
encing the work of Kac and van Moerbeke [13], Flaschka [9] and Moser [15], [16], 
see also Moser and Zehnder [17]. The technique follows a change of variables 

(3.17) Uk = e'" 

under which the differential equation (|3.ip takes the form 

(3.18) ^ ^ (e^-^-i _ e^-'+M, 

da 

which we relate to the equation 

(3.19) ^ = (e^'=-i-^'= + gXfc-x.+i) 
via the relation 

(3.20) Zk=Xk-Xk+i. 

Observe that due to the periodicity of the x^s the sum of the z^s is zero, which 
might be a source for an extra invariant. Equation p.l9p is the well studied Toda 
lattice that forms a completely integrable Hamiltonian system. It is clear, however, 
from (|3.19p . that the solutions we are interested in, namely, solutions of (|3.18p . 
are differences of unbounded monotone increasing solutions (Moser [16] computes 
the asymptotics of these solutions). Due to the unboundedness, this Hamiltonian 
structure is not of much help in the computations of solutions of (|3.ip . We do not 
know if, possibly, another change of variables will induce a Hamiltonian structure 
on p.ip or will make solutions of (|3.ip images of bounded solutions of a Hamiltonian 
system. 

4. A COMMENT ON THE FAST-SLOW SYSTEM 

In this section we show that the positive orthant is invariant with respect to the 
full fast-slow system (11.11) . Moreover, we show that every such solution converges 
to a constant state Ui ~ U2 = ■ ■ ■ = Um- 

Proposition 4.1. System hl.l]) with Uk{0) > for all k = 1,2, . . . , N , possesses 
a unique positive solution Uk{t) > 0, for all t > 0. Moreover, 

N N 

(4.1) l[Uk{t)>l[Uk{0), 

k=l k=l 



8 Z. ARTSTEIN, C.W. GEAR, I.G. KEVREKIDIS, M. SLEMROD, AND E.S. TITI 



and Uk(t) converges to a constant, independent of k, for all k — 1,2, .. . , iV, as t 
goes to infinity. 

Proof. Since Uk{0) > 0, for all fc = 1. 2, . . . , A^, then for a short time we conclude 
from (fLTj) that 

-^logt/.w =^E[^4 

k=l k=l 
A.;— 1 

> 0. 

The last inequality is true since the function f{U) = is monotone decreasing. 
Hence, 

AT JV 

(4.3) ^logC/fe(t) >^logC/fe(0), 

fc=l A:=l 

and 

JV JV 

(4.4) llUk{t)>l[Uk{0)>0. 

fe=i fe=i 

Consequently, C/fc(i) > for all fc = 1, 2, . . . , N, and for all < > 0. 
Since in addition 

JV JV 

(4.5) 5]C/fe(t) = ^C/..(0), 

k=l k=l 

we have Uk{t) uniformly bounded and hence we have global unique solutions. Fi- 
nally, the expression — E^j^ log Uk{t) is a Lyapunov function and via the well known 
LaSalle Invariance Principle (see, for example, [12], Chapter 10, Theorem 1.3, page 
316) all solutions converge to the largest invariant set contained in the set defined 
by the solution of the algebraic equation 

^ 1 1 

Since the only possible solutions to the above equality are Ui — U2 = ■ ■ ■ = Un = 
const, the proposition is proven. 



5. Young measures solutions and observables of the fast-slow 

PROBLEM 



In this section we address the full system p.7p . Let T > be arbitrary, but fixed. 
Based on the theory developed in [3] we represent Hmits, as e ^ 0, of solutions as 
Young measures over the interval [0, T]; then examine how slow observables of the 
system evolve over this arbitrary, but fixed, interval of time. The analysis is the 
basis for the computational developments in the closing section, section 6. 

In the discussion that follows N is fixed, hence h is fixed. We will write U for 
the vector (f/i, U2, ■ ■ ■ , Uk) in . We are interested then in limits of solutions 
of (|27l) . as e ^ 0, over the interval [0, T]. 
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The Young measures we use are defined on a time interval [0,T], with values 
being probability measures on , namely, on the state space of solutions of (I2.7|l . 
For the relevant theory of Young measures see, e.g., [3] and references therein. A 
solution of (|2.7|) . say Us{t) : [0,T] , can be viewed as a Young measure 

when for each r the vector U^{t) is interpreted as a Dirac measure supported 
on {C/e(r)}. A Young measure can be viewed also as a measure on [0,T] x R^ . 
Convergence of Young measures of the full fast-slow system (|2.7p . as e — > 0, is 
interpreted as the weak convergence of the associated measures on [0, T] x i?^; the 
notion includes the case where a sequence of point-maps into i?^, interpreted as 
Young measures, converge to a Young measure. Convergence of solutions of (|2.7p 
in the sense of Young measures reflects that the distribution of the location of the 
solutions converge. In particular, the convergence in the sense of Young measures 
gives a more complete description of the limit behavior of the full system than 
the standard averaging method. This is because the Young measures approach 
provides a mathematical framework, i.e. in the sense of measures, for investigating 
and describing the limit behavior of both the slow as well as fast dynamics, while 
the averaging method restricts itself to the slow part of the dynamics. 

The following result addresses the convergence of solutions in the sense of Young 
measures; its proof is based on the developments in previous sections and on [3], 
Theorem 4.4. 



Theorem 5.1. Let the initial condition for \2. 7\ ) be given (alternatively, let 

be in a compact set in the positive orthant of W^). Let T > be given, and 
fixed. Denote by Us{t) the solution of {2. 7[ ), for a given e, over the interval [0,T]. 
Then a subsequence — > exists such that U^^ (t) converge as k ^ oo, in the sense 
of Young measures on [0,T], to a Young measure, say Hq{t), t G [0,T]. Moreover, 
for almost every t e [0,T] the measure fJ-oir) is an invariant measure of the fast 
equation i3.1]} . 

Proof. A straightforward computation shows that the quantity Ui + ■ ■ ■ + Un is 
not only an invariant of the fast flow (|3.ip , see (|3.12p , but also of the full equation 
(2.7) (indeed, it represents the total mass of the system, which, under periodic 
boundary conditions, is preserved under viscosity). Hence the family Ueir) of 
functions, which is from [0,T] into the positive orthant of i?^, is bounded in R^ . 
In particular, each solution in the sequence can be extended to the entire interval 
[0,T]. Now the theorem follows from [3], Theorem 4.4. 

The next result addresses the evolution of observables of the system. A measure- 
ment, or observable, is simply a function of U. The notion of an orthogonal slow 
observable was deflned in [3]; in our system it amounts to a measurement which 
is a flrst integral, namely, an invariant, of the fast equation p.ip . In Section 3 we 
identifled + ^ such invariants. Denote by Vj{U),j = 1,2, . . . , ^ the invariants 
which are the ^ traces of LP,p — 2,4, . . . , N (see Section 3, but note that in that 

section the traces are expressed in terms of Ak = U^); denote by vn_^i{U) the 
invariant given in ()3.15p . The invariance of each of the Vj implies that whenever fi 
is an invariant measure of (|3.ip . the observed value Vj(U) is constant on the support 
of /i. We denote this value by VjifJ.)- 

Theorem 5.2. Let U^^ (r) be as in the statement of Theorem 5.1 which converge, 
as k oo, in the Young measures sense, to the Young measure IJ-o{t), for r £ 
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[0,T]. Denote Vji^r) — i)j(/io(T)), namely, the measurement on the limit invariant 
measures. Then for each Vj,j = 1, 2, . . . , + 1, the function 'Oj(/io(T)) satisfies the 
differential equation 



where G{U) = G{Ui, . . . , Un) is the vector field given in i2.8\) . and the V operator is 
with respect to the vector U . Furthermore, the sequence of measurements Vj(Ue^ (t)) 
converge weakly to Vj{iio{T)), as fc — > oo, for j = 1, 2, . . . , + 1- 
Proof. The claims form a particular case of Theorem 6.5 in [3]. 

Under stronger conditions on continuity properties of the limit measure fJ-oir), 
the weak convergence in the previous result becomes a strong convergence; we do 
not pursue this here since in our system it is not easy to establish such continuity. 

The differential equation (|5.ip is not an autonomous ordinary differential equa- 
tion since the limit measure fJ-oir) is not determined by the observation Vjir). In 
some cases a set of observables determines the invariant measure; then, under con- 
tinuity assumptions one can write an ordinary differential equation for this set of 
observables; without continuity a differential inclusion would determine the limit 
evolution of the observables. See Theorem 6.9 and Remark 6.11 in [3]. These re- 
sults can be stated in our case employing the terms in (|5.ip , (12. 8p ; we do not follow 
this direction since we do not know to what extent the ^ ^ ^ observables we found 
form such a determining vector of observables. At any rate, the structure described 
here allows one to make progress toward a computational method, as displayed in 
the next section. 



In section we present a computational study of the theoretical tools that were 
developed in [3] and which were presented in the previous sections for system p.7p 
(or equivalently (|2.8[) ). with small values of e. Specifically, we introduce two ef- 
ficient numerical schemes for computing the evolution of the slow observables of 
(|2.7p . for small values of e. In particular, we implement these schemes for com- 
puting the evolution of the slow observables Vj{U), j ~ 1, 2, . . . , y + 1, that are 
discussed in section 5. Moreover, we demonstrate the computational efficiency of 
these algorithms in comparison to direct numerical simulations of the full system 
p.7p for small values of e. Both schemes are based on numerical evaluation of the 
time derivatives of the slow observables. That is, they are based on evaluating the 
time derivatives of the invariants of the fast system (|3.ip . Vj{U), j = 1, 2, ■ • ■ , ^ + 1' 
as they are drifted by the slow field G{U) in (|2.8p . 

In addition to the evolution of the slow observables, the method can simulate 
the evolution of the invariant measures of (|3.ip themselves, thus getting a picture 
of the entire rapidly oscillating trajectories of (|2.7p . 

The first approach, the Young measures approach, is based on computing the 
right hand side of (|5.ip . To this end, one has to compute first the invariant measures 
HoiT) (see Theorem 5.1 and Theorem 5.2) of the fast system p.ip by, for example, 
integrating system p.ip for an appropriate interval of time, and then to advance 
Vj{T), j = 1, 2, . . . , Y + Ij numerically, according to the numerical evaluation of the 
right hand side of (|5.ip . 



(5.1) 




6. Computational study 
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The second approach, the equation-free approach of [14] (it is equation-free since 
it avoids the direct use of any analytic form for the equation of evolution of the slow 
observables Vj), is based on numerically evaluating the time derivative of Uj(L/j(r)), 
j = 1, 2, . . . , Y + I7 from the detailed simulator of the solution Ueir) of the full 
fast-slow system (I2.7p . for small values of e. It is worth mentioning, however, that 
the success of this method relies on the a priori knowledge that Vj{Ue{T)), are slow 
observables of (|2.7p and that at the limit, when e ^ 0, their evolution is governed 
by a deterministic equation; that may not be at hand. These facts were established 
in the previous sections. 

The efficiency of both algorithms stems from the fact one can take relatively 
large time steps for advancing the slow observables, Vjir) in the first approach, and 
Vj{Ue{T)), for small values of e, in the second approach. We call this step in both 
algorithms the projective step. 

It is worth mentioning the philosophical difference between the two approaches. 
The first approach, which is based on Theorem 5.1 and Theorem 5.2, takes ad- 
vantage of the fact that the evolution of the slow observables ^^([/^(r)), j = 
1,2,...,^ + 1, along the solution [/^(t) of system ()2.7p . are approximated, for 
small values of e, by their corresponding limits %(t), j = 1, 2, . . . , ^ + 1j as e ^ 0. 
Therefore, the first approach focuses on computing the slow observables 'Oj(t) of 
the limit equation, as e ^ 0. This is because the slow observables Vjir) will then 
be considered as the approximations to the corresponding "real" slow observables 
Vj{Ue{T)), for small values of e. The second approach, which is based on [14], 
deals, on the other hand, directly with system (|2.7p without resorting to the limit 
equation, equation (jS.ip . but it implicitly assumes the existence of a limit equation. 

Recall that in section 3 we established that there are at least + 1 invariants 
of the fast system (|3.1[) . therefore each trajectory is constrained to a manifold of 
dimension at most y — 1- This appears to be achieved for the case of TV = 6 
as shown in Figure [Tl which illustrates the solution of (|3.ip for one initial value 



3^^ 



2.5 




Figure 1 . Torus for the case = 6 of system p.ip . Initial values 
were [1113 2 1]. 

plotted in the first three variables, Ui{a),U2{cr) and U3{a). The solution hcs on 



12 Z. ARTSTEIN, C.W. GEAR, I.G. KEVREKIDIS, M. SLEMROD, AND E.S. TITI 



the surface of a 2-torus in this 3D projection. The same observation holds for all 3D 
projections of this case except those projections into either of all the odd-numbered 
or all the even-numbered coordinates (where it necessarily lies on the surface of 
one sheet of a two sheeted hyperboloid; this is because of the invariants noted in 
Observation 3.1). The trajectory "orbits" around the "long" direction of the torus, 
precessing from orbit to orbit. We assume, however, that we can approximate it, 
and the relevant invariant Young measure, /xq, by integrating over an appropriate 
finite interval of time. Indeed, we have found that integrating for a few "orbits" 
over the "long" direction of the torus can give an adequate approximation for our 
application. 

For = 6 a typical solution of ()3.1|) was found to be approximately a combina- 
tion of two periodic motions, the fastest corresponding to an orbit around the long 
direction of the torus and the slow one corresponding to the precession time around 
the short direction. See, for instance, the graph of the first component, Ui{(j), of 
the solution of system (|3.ip shown in Figure [2l 



3.5 




Figure 2. Ui{a) for the case in Figure [T] 



Whereas the general theory, as described in section 3, shows that there is a set 
of at least ^ + 1 invariants for all solutions of the iV-dimensional problem (|3.ip . 
some particular solutions may satisfy additional invariants which we might call local 
invariants. That there are initial conditions leading to local invariants is easy to 
see. For instance, the case N — 3 has 2 invariants {Ui + U2 + and U1U2U3). 
The case = 6, the solutions of (|3.1[) that have started with Ui = U4, U2 = C/5, 
and C/3 = Uq, will continue to maintain those identities along the trajectory and 
will have the same solution as the = 3 problem with those starting conditions. 
Hence its solution will lie on a one-dimensional manifold. This is shown in Figure 
131 (The invariants in this case are the 2 invariants, for the = 3 case, and the 3 
equalities, Ui — Ui+3.) The case A^ — 12, for example, that is known to have at least 
7 invariants, will have certain initial conditions with 10 invariants (corresponding 
to 2 copies of the A^ = 6 case) or 11 invariants (corresponding to 3 copies of the 
A^ = 4 case, or 4 copies of the A^ = 3 case). 
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Figure 3. Case iV = 6 as pair of = 3 cases for system p.ip . 
Initial values were [3 2 1 3 2 1]. 



6.1. Two approaches for computing slow observables. In this subsection we 
illustrate in more detail the two approaches for computing the behavior of the 
fast-slow system (j2.7p which we rewrite in the form 

(6.1) ^ = f/fe(f/fc-i - C/fc+i) + iy{Uk+i - 2Uk + Uk-i), for fc = 1, 2, • • • , TV, 

where ^ — very small. 

We do this for the case iV — 6, but the technique is not restricted to that 
case. Both approaches are based on evaluating the time derivatives of the slow 
observables, Vj{U), j = 1,2, ■ ■ ■ , ^ + Ij as we have discussed above. 

In the first approach, which is based on the Young measures as described in 
section 5, we utilize the right hand side of equation (|5.ip to evaluate the "slow" 
derivatives of the slow observables f)j(r), j = 1, 2, • • ■ , + 1. This can be done 
by integrating the fast system (|3.ip . namely, the v — case of equation (|6.ip . to 
approximate the invariant Young measure /io of p.ip : what we practically do to 
evaluate the right hand side of (|5.ip is to integrate over equally spaced time steps 
and approximate the integral in the right hand side of (|5.ip by an averaged Riemann 
sum based on those points. 

In the second approach, the equation-free method [14], we integrate the actual 
full fast-slow system equation (|6.ip with small value of > 0, compute the slow 
observables Vj{U„{t)), j = 1, 2, • • • , ^ -I- 1, along the trajectory U^{t), which are 
guaranteed to be slow by our theory, and evaluate their slow derivatives by numer- 
ical differencing. 

In both approaches we use the time derivative evaluation of the slow observables 
Vj , j = 1,2, ■ ■ ■ , Y + 1, in order to advance them in time. We call this very last step 
in our scheme the "projective level" . We observe that the gain in the speedup of 
both schemes is due to the use of relatively large time steps at the projective level 
of the schemes for advancing the slow obervables. In order to repeat the process 
we need to initialize the detailed simulator (i.e., find initial conditions, U, for the 
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variables U in the original iV-dimensional space) consistent with the "projected" 
values of the ^ ~^ ^ observables. That is, from the computed values of the 
slow observables, say vj, = 1, 2, • • ■ , + 1, we need to find a corresponding state 
variable U S i?^, for which Vj{U) = vj, for all j = 1, 2, • • • , + 1. This procedure 
is called "lifting" in the equation-free literature [14]. In general this will involve 
the solution of a set of ^ + 1 nonlinear algebraic equations with N unknowns 
(and this implies that — 1 features of the solution can, in principle, be selected 
(almost) arbitrarily; in our six-dimensional example this "arbitrary" selection was 
implemented by carrying the values of two of the entries of the sought vector U 
over from the last simulation.) As mentioned above, we take the point U € i?^ 
as an initial value for the corresponding system, i.e., an initial value for system 
(|3.1[) to compute the corresponding Young measure in the first approach, or an 
initial value for (|6.ip to compute the corresponding solution Ui,{t) in the second 
approach, and we repeat the above described process. Note that if t/ G is 
another point that has the same projected values of the slow observables as those 
of U, i.e. Vj{U) = Vj{U) = Vj, for all j = 1, 2, • • ■ , -1-1, then we are not guaranteed 
that the dynamics of the fast system ()3.ip starting from U will lead to the same 
invariant Young measure, /zq, which is corresponding the trajectory starting from 
U (see the relevant discussion regarding this matter in the closing paragraph of 
section 5). This is unless the invariant Young measure of (j3.1l) is unique, which is 
the case if the trajectory corresponding to the solution of fast dynamics p.ip "fills" 
the manifold defined by the invariants of (j3.ip . As noted earlier, this manifold can 
be up to — 1 dimensional, since it is a sub-manifold of the manifold defined by 
the intersection of the invariant manifolds Vj{U) = constant, j = 1, 2, • • • , ^ -I- 1; 
although in reality it might be of lower dimension because of the additional local 
invariants. 

6.2. The first approach - using the invariant Young measures. The invari- 
ants Vj(U), j = 1, 2, • ■ • , Y + of the fast system p.ip were identified in section 3 to 
be slow observables of the full fast-slow system (|2.7p , or equivalently of (|2.8p . Since 
these quantities are constants of motion for the fast dynamics, the contribution of 
the fast vector field, £~^F{U) in (|2.8p . to their evolution is zero. Consequently, we 
focus only on computing the evolution of these slow observables as they are drifted 
by the slow vector field G{U) in (|2.8p . At the limit case, as e ^ 0, the rate of 
change of these slow observables, -§^Vj, is given by the right hand side of (|5.ip . that 
is, by the average of Wvj{U) ■ G{U), the directional derivative of Vj{U) along the 
vector field G{U), with respect of the Young invariant measure of the fast dynamics 
p.ip . The need to average Vvj{U) ■ G{U) over the measure can be seen from the 
wild fluctuation in Figure |4]of Vvz{U{(j)) ■ G{U{cr)), for the 3rd Lax invariant V3, 
along the trajectory U{a) of the fast dynamics (|3.ip that is shown in Figure [TJ 

As we mentioned earlier in the introduction of section 6, a typical solution U{cr) of 
the fast dynamics p.ip was found to be approximately a combination of two periodic 
motions, one fast and the other slower. The average value of Vw3([/(cr)) ■ G{U{a)) 
along the trajectory U{a) of the fast dynamics (|3.ip . from the time a = to the 
current time, is shown in Figure [5] Either the time over which the average is taken 
has to be chosen carefully with respect to the approximate period of U{a), or the 
average must be taken over a long time. In Figure[5]the circles indicate the averages 
at "periodic" points (the approximate fast period was estimated as 2.4868). Those 
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Figure 4. Evolution of Vu3(C/((t)) ■ G{U{a)) along the trajectory 
U{a) of Figure [H 

dots indicate that we can get a reasonable approximation for the Young measure 
using a small number of approximate orbits. 
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Figure 5. Averaged Vv3{U{a)) ■ G{U{a)) along the trajectory 
U{a) of Figure [H 

We use this method of derivative evaluation of the slow observables, by virtue 
of the average over the Young's measure in the right hand side of (jS.ip . in order 
to integrate them forward in time by simple forward Euler integration scheme with 
relatively large time steps. Then we compare the performance, i.e. speedup and 
efficiency, of this method to computing the slow observables through direct numer- 
ical simulation of the full fast-slow system (|6.ip with i/ = 10"'*, where the latter 
will be considered as a "true" solution. Computing the slow observables via the 
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Young's measure approach was done separately by first evaluating the time deriva- 
tives of the slow observables, as it was described above, over a fixed number of "fast 
periods" , and then advancing the slow observables by taking an integration step in 
a forward Euler scheme of either 3, 6, or 12 period length. We emphasize that the 
gain in the efhciency of this method lies in the last step, i.e. in the use of these 
relatively large time steps in advancing the slow observables forward in time. The 
solutions for the second through fourth slow observables, V2,V3,V4, over 3,000 "fast 
periods" (a time interval of 7,460.4) are shown in FigurelHl The slow observables of 
the "true" solution Ui,{t) (found by integrating the full fast-slow system (|6.ip . with 
f = 10"'', and then computing the slow observables Vj{U^{t)) directly from it) are 
plotted as solid lines, while the integration of the slow observables via the Young's 
measure method, as described above, are plotted as dotted lines; but the difference 
(to plot accuracy) is not visible. We ran separate integrations over either 1, 2, 
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Figure 6. Behavior of the slow observables V2,v^ and as they 
are drifted by the slow diffusion. 

3, or 4 fast periods, in order to approximate the average in the right hand side 
of equation (jS.ip with respect to the corresponding approximate invariant Young 
measures of the fast system (|3.ip . These separate integrations, or approximations 
of the Young measures, had less than a 10% effect on the solution errors. The 
dominant errors were due to the Euler integration. The errors in the slow observ- 
ables after 600 fast periods are shown in Table [T] The speedup and efficiency of the 



Euler step (periods) 


V2 


W3 




3 


0.45 


0.47 


0.68 


6 


1.3 


2.1 


1.8 


12 


3.9 


10.8 


4.1 



Table 1. Errors (scaled up by 10^) at t = 1,492.08 (600 fast 
periods) in the slow observables integrated using approximations 
to Young measure averages. 
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Young measure integration scheme, in comparison to direct numerical simulation 
of the full fast-slow system ()6.ip . depends on the value of v. The cost of evaluating 
the time derivative of the slow observables via the Young measure averaging in the 
right hand side of (|5.ip is determined solely by the accuracy needed. Moreover, 
the number of integration steps needed for advancing the slow observables forward 
in time is also determined solely by the accuracy needed. However, if we were to 
integrate equation (|6.ip directly, the number of time steps needed would be propor- 
tional to the integration interval which is proportional to l/v. If, for example, we 
needed 50 integration steps for one approximate period, v = 0.001, and we wanted 
to track the dissipative (diffusion) terms until they had reduced by a factor of 10"^ 
we would need to integrate from t = {) to t = 11,500 (approximately). Since each 
approximate period is 2.4868 we would need about 4,600 periods, or 23,000 inte- 
gration steps for a direct numerical integration of ()6.ip . If we used a higher-order 
integration method (rather than the first order Euler method used in this illus- 
tration) at the projective level, namely for advancing the slow observables in time 
after evaluating their time derivatives, we could probably integrate the decaying 
slow observables with much less than 100 projective steps - each of which requires 
one integration around of the fast system using 50 steps, for a total of 5,000 steps, 
or a speedup of better than 4 to 1. If were reduced by one tenth, the speedup 
would increase to better than 40 to one. 

As mentioned, while computing the invariant measures of the fast system p.ip 
(or equivalently of (j6.ip with v ^ Q) we actually simulate the entire evolution of 
the limit system of ()6.ip . as v ^ 0, via the supports of these measures - tori in our 
example. The decay of the tori of the full fast-slow system (16. ip . with > 0, can 
be seen in Figure [7] which shows the trajectories after 600, 1200, 1800, 2,400, and 
3,000 fast periods of integration. 




Figure 7. Evolution of Tori of system (|6.ip for initial condition 
[11114 1]. 
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6.3. The second approach - the Equation-Free method. The equation-free 
method [14] utilizes the explicit knowledge of the numerical solution U,^{t) of (16. 1|) . 
and consequently of the observables Vj{U,^{t)), j = 1, 2, • • • , ^ + Ij for short interval 
of time in order to approximate, here by numerical differencing, the time deriva- 
tive of Vj{Uiy{t)) in order to project forward their values. The method is called 
equation-free since it does not use explicitly any equation of motion for the slow 
observables. However, it assumes that the evolution of Vj{Uu{t)) is slow, and that 
there is an underlying process, or equation, that governs their evolution, which 
we do not use explicitly. (A valuable research direction would be to develop a 
method that identifies, possibly numerically, the slow observables in case they are 
not prescribed.) 

For small values of i/ > the evolution of some of the slow observables, Vj{Uu{t)), 
j = 1,2, ••■,^-1-1, along solutions U^{t) of (|6.ip can exhibit a slight oscillatory 
behavior. For example, Figure [5] shows the early behavior of v^{U^{t)) along the 
solution U^{t) of equation (|6.ip for the initial condition [11114 1]. It exhibits a 
fine-scale oscillation. However, the values at the "periodic" points are reasonably 
smooth so we can use those points in order to evaluate the time derivative of the 
slow observables, VjiUuit)), by numerical differencing. Then we use the numerical 
time derivatives of Vj{U^(t)) to project them forward in time. The results are 
illustrated in Figure [9l We show the results of projective integration of the slow 
observables (dotted lines) versus the evolution of the slow observables along a true 
solution of (|6.ip (solid line) for the three slow observables U2, and v^. Again, and 
as in the Young measure method, the results are indistinguishable to plot accuracy. 
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Figure 8. Slow observable V3{U„{t)) for initial condition [11114 1]. 



Table [2] shows the errors in the slow observables at t = 1, 492.08 using projective 
Euler integration by evaluating the time derivative from the slope of the chord over 
one period. 
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Figure 9. Behavior of the slow observables V2^vz and V4 as they 
are drifted by the slow diffusion using the Equation-free integra- 
tion. 



Euler step (periods) 


V2 


V3 


V4 


3 


0.8 


3.6 


0.5 


6 


1.9 


5.9 


1.8 


12 


4.2 


14.2 


3.9 



Table 2. Errors (scaled up by 10^) at t = 1,492.08 (600 fast 
periods) in the slow observables integrated using the Equation-free 
method over one period. 
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